/*
 * EllipsePore_test.cpp
 *
 *  Created on: Nov 14, 2011
 *      Author: yi
 */

#include "mbuild/EllipsePore.h"
#include "mbuild/DxFile.h"
#include "mbuild/Tfunc3D.h"
#include "time.h"
#include "mbuild/DefaultFunctions.h"
#include "mbuild/Tfunc3D.h"
#include <iostream>

int testSetupPore(EllipsePore& pore) {
	double R_in,R_out,T;
	double xmin,xmax,ymin,ymax,zmin,zmax;
	R_in=10;R_out=20;T=27;
//	cout<<"R_in:";
//	cin>>R_in;
//	cout<<"R_out:";
//	cin>>R_out;
//	cout<<"Membrane Thickness(T):";
//	cin>>T;
	try {
		pore.setup(R_in,R_out,T);
	} catch(ParameterError& error) {
      cout<<"Incorrect Parameters!\n";
      return 0;
	}
	xmin=ymin=zmin=-40;xmax=ymax=zmax=40;

//	cout<<"Xmin:";
//	cin>>xmin;
//	cout<<"Xmax:";
//	cin>>xmax;
//	cout<<"Ymin:";
//	cin>>ymin;
//	cout<<"Ymax:";
//	cin>>ymax;
//	cout<<"Zmin:";
//	cin>>zmin;
//	cout<<"Zmax:";
//	cin>>zmax;
	//Time to construct distance map. Reducing this: calculate only the pore region?
	//More accuracy? larger R in fast matching? using df/dx,df/dy,df/dz?
	clock_t t1=clock();
	//TODO: Problem: xmin... has to be changed if region > box size to increase efficiency.
    pore.set_region(xmin,xmax,ymin,ymax,zmin,zmax,true,0.5,0.5,0.5);
	clock_t t2=clock();
	cout<<"Finished constructing distance map in:"<<double(t2-t1)/CLOCKS_PER_SEC<<" s.\n";
    double cutoff;
    int ncyc;
    cutoff=10;ncyc=3;
    //TODO : setup calculation.

//    cout<<"Integration cutoff(for improving the accuracy of calculating density around the surface):";
//    cin>>cutoff;
//    cout<<"Integration cycles(more cycles,higher accuracy but longer calculation time):";
//    cin>>ncyc;
    try {
    pore.set_cutoff(cutoff);
    pore.set_integration_cyc(ncyc);
    } catch(ParameterError& error) {
    	return 0;
    }
    double r0=11,r1=11;
    //TODO : add cylindrical pore.
//    cout<<"Has cylindrical pore?\n";
//    cout<<"R up:";
//    cin>>r0;
//    cout<<"R bottom:";
//    cin>>r1;
    CylindricalHole hole=CylindricalHole(r0,r1,T+10*cutoff,-T/2-5*cutoff);
    if (r0 >0 && r1 >0) {
    pore.has_hole(&hole);
    }
	return 1;
}

int testConstruct(GridSpace& space,EllipsePore& pore) {
	Gaussian g;
	g.set_offset(1.0);
	g.set_width(0.75);
	g.set_sigma(1);
	g.Solute_Do_Not_Eql(0.0);
	clock_t t1=clock();
	pore.set_distribution(&g);
	GridSpace ref;
	space.copyto(&ref);
	pore.fill(&space,&ref);
	clock_t t2=clock();
	cout<<"Finished filling Gaussian in:"<<double(t2-t1)/CLOCKS_PER_SEC<<" s.\n";
	pore.create_hole(&space,&ref);
	clock_t t3=clock();
	cout<<"Finished Creating hole in:"<<double(t3-t2)/CLOCKS_PER_SEC<<" s.\n";
	DxFile::write_dx("add_charge.dx",&space);
	return 1;
}
int testConstructSwitch(GridSpace& space,EllipsePore& pore) {
	Switch s;
	s.set_values(0.0,10.0,80.0);
	s.set_width(5.0);
	s.Solute_Eql(2.0);
	clock_t t1=clock();
	pore.set_distribution(&s);
	GridSpace ref;
	space.copyto(&ref);
	pore.fill(&space,&ref);
	clock_t t2=clock();
	cout<<"Finished filling Switch in:"<<double(t2-t1)/CLOCKS_PER_SEC<<" s.\n";
	pore.create_hole(&space,&ref);
	clock_t t3=clock();
	cout<<"Finished Creating hole in:"<<double(t3-t2)/CLOCKS_PER_SEC<<" s.\n";
	DxFile::write_dx("add_diel.dx",&space);
	return 1;
}
int testConstructGouy(GridSpace& space,EllipsePore& pore) {
	Gouy g;
	clock_t t1=clock();
	g.set_lipid(65,1,1); //Area per lipid 65A^2,anionic lipid fraction 0.3, net charge per lipid 1.
	g.set_ion(1,0.1);
	g.set_offset(3);
	g.set_temp(298);
	pore.set_distribution(&g);
	GridSpace ref;
	space.copyto(&ref);
	pore.fill(&space,&ref);
	clock_t t2=clock();
	cout<<"Finished filling Gouy in:"<<double(t2-t1)/CLOCKS_PER_SEC<<" s.\n";
	DxFile::write_dx("gouy.dx",&space);
	return 1;
}


int testInhomogenuousConstruct(GridSpace& space,EllipsePore& pore) {
	Gaussian g;
	g.set_offset(3.0);
	g.set_width(5);
	g.set_sigma(1);
	//TODO add inhomogenous function
//	Tfunc3D f;
	pore.set_distribution(&g);
	GridSpace ref;
	space.copyto(&ref);
	pore.fill(&space,&ref);
	pore.create_hole(&space,&ref);
	return 1;
}
int main(){
	EllipsePore pore;
	testSetupPore(pore);
	GridSpace space;
	DxFile::read_dx("charge_X.dx",&space);
	testConstruct(space,pore);
	DxFile::read_dx("dielx_X.dx",&space);
	testConstructSwitch(space,pore);
	space.clear(0.0);
	testConstructGouy(space,pore);
	return 0;
}
